Self-sustained spatiotemporal oscillations induced by membrane-bulk coupling 
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We propose a novel mechanism leading to spatiotemporal oscillations in extended systems that 
does not rely on local bulk instabilities. Instead, oscillations arise from the interaction of two sub- 
systems of different spatial dimensionality. Specifically, we show that coupling a passive diffusive 
bulk of dimension d with an excitable membrane of dimension d — 1 produces a self-sustained os- 
cillatory behavior. An analytical explanation of the phenomenon is provided for d = 1. Moreover, 
in-phase and anti-phase synchronization of oscillations are found numerically in one and two dimen- 
sions. This novel dynamic instability could be used by biological systems such as cells, where the 
dynamics on the cellular membrane is necessarily different from that of the cytoplasmic bulk. 

PACS numbers: 82.40.Ck, 87.16.Ac, 47.54.+r 
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Introduction. In spatially extended systems, spa- 
tiotemporal oscillations usually arise from short- 
wavelength, finite-frequency instabilities that affect the 
local dynamics of the system's bulk. Within that sce- 
nario, boundaries are reduced to passive elements that 
play somewhat secondary roles, such as wavelength dis- 
cretization and wavevector selection [l| . There are many 
situations in nature, however, where boundaries have ac- 
tive dynamics. Fronts are known to be initiated, for in- 
stance, at the interface between different catalytic com- 
ponents in microcomposite surfaces [2[ . Similarly, chaotic 
dynamics has been shown to arise in a catalytic surface 
coupled to a (passive) gas phase Other examples 
of active surfaces include Langmuir monolayers Q and 
membranes with active proteins such as proton pumps 
Q . These systems delimit regions of higher dimension- 
ality, which usually have different dynamics from that of 
the active boundary. 

Few studies have addressed the interplay between dif- 
ferent dynamics of a bulk and its boundary. In [^] it 
has shown that an active membrane can give rise to 
stationary cytoplasmic patterns. Special attention has 
been paid to the issue of pole-to-pole protein oscillations 
underlying symmetrical cell division in bacteria 0, 
Most models of this phenomenon assume that the main 
source of the oscillations are biochemical reactions oc- 
curring at the cell membrane, with the cytoplasm being 
mainly passive, hosting at most phosphorylation reac- 
tions d Eil El, e m. It is thus of our interest to 

determine whether nontrivial dynamics can arise in the 
limiting case of an active boundary delimiting a purely 
passive bulk. 

Motivated by this system, which is known to include 
activator and inhibitor proteins [7j, and taking into ac- 
count the fact that activator-inhibitor dynamics sustains 
excitability 1J], we will consider in what follows an ex- 
citable membrane, limiting an otherwise purely passive 
bulk [l5j |. Indeed, we will show that this simplified sce- 
nario is able to sustain dynamic spatiotemporal oscilla- 



tions in a wide parameter range, even though neither 
the bulk nor the boundary are oscillatory. Moreover, we 
will provide an analytical explanation for this effect in 
the case of a one-dimensional bulk and point-like bound- 
aries. This analysis allows us to predict the extent of the 
oscillatory region in terms of the relevant parameters: 
the system length and the coupling strength between the 
bulk and the boundaries. 

The model. We consider a spatially-extended passive 
system affected by simple diffusion and linear degrada- 
tion, bounded by an active membrane with activator- 
inhibitor dynamics. The equations that mathematically 
describe the bulk are 



d t U = DuV 2 U - auU, 
d t V = D V V 2 V - a v V, 



(1) 

(2) 



where U and V are the concentrations of activator and 
inhibitor, respectively, the corresponding diffusion coeffi- 
cients are Djj and Dy, and both species are assumed to 
decay at rates ajj and ay [3]. The results shown below 
do not change in the presence of a constitutive expres- 
sion of the bulk species, represented by the addition of 
constant terms in the r.h.s. of equations (HJ)- ([2J) - 
The dynamics at the system's boundary is given by 



u = f(u, v) + k u (n- VU), 
v = eg(u, v) + k v (h ■ VV). 



(3) 
(4) 



The reaction terms f(u,v) and g(u,v) are chosen to ac- 
count for a local activator-inhibitor dynamics, so that the 
w-nulcline [f(u,v) = 0] has a cubic shape in the (u,v) 
space, while the w-nulcline [g(u,v) = 0] is monotonically 
increasing, as in typical Fitzhugh-Nagumo models (see 
Fig. [2] below). In particular, we have considered the fol- 
lowing expression in dimensionlcss form 



f(u,v) =u-q{u 
g(u, v) = uz — v, 



(5) 
(6) 
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where z and q control the shape of the nulclines. The pa- 
rameter e that accompanies g(u, v) determines a different 
time scale for both species, so that u is much faster than v 
if e -c 1. We consider in what follows z = 3.5 and q = 5, 
which renders the membrane excitable. Under this con- 
dition, the membrane (when isolated from the bulk) is 
at rest, and only when a small perturbation is applied 
to it, an excursion corresponding to an activator pulse 
is produced. The second term in the right-hand-side of 
Eqs. ([3|) and (jlj accounts for the exchange of species be- 
tween the membrane and the bulk. The constants k u and 
k v determine the coupling strength, while n-VC(r) is the 
normal derivative of the concentration field C (r) , with h 
being a unit vector normal to the boundary and pointing 
towards the bulk. 

Phenomenology. In principle, one would expect the 
system described above to be quiescent unless a per- 
turbation is applied to the membrane. Such a pertur- 
bation would excite a concentration pulse at the mem- 
brane, which would then propagate into the diffusive bulk 
and progressively decay. However, coupling between the 
membrane and the bulk can give rise to interesting new 
phenomenology. Due to the coupling, degradation of the 
inhibitor in the bulk leads to a decrease in the concen- 
tration of the inhibitor also in the membrane, which sub- 
sequently allows, via the activator-inhibitor dynamics, a 
pulse in the activator concentration. After some time, 
the excitability mechanism reduces the activator level 
spontaneously back to the resting state, and the process 
can restart again and repeat endlessly, leading to self- 
sustained oscillations even when neither the membrane 
nor (evidently) the bulk are intrinsically oscillatory. This 
is indeed observed in our model, as we show below. 

Once oscillations are autonomously occurring, we can 
analyze what happens when the bulk is limited by two 
opposing boundaries, similarly to an ellipsoidal bacterial 
cell. If the distance between the poles is small enough 
the oscillations interfere, and eventually become synchro- 
nized. This is shown in the one-dimensional numerical 
simulations presented in Fig.[T] Two regimes correspond- 
ing to in-phase and anti-phase oscillations of the two 
poles are observed. For each regime, a plot of the time 
evolution of v at the two boundaries is displayed. Addi- 
tionally, spatiotemporal plots are shown in the right pan- 
els. The results clearly show not only that self-sustained 
oscillations appear even when the system is not intrinsi- 
cally oscillatory, but that the oscillations can synchronize 
either in phase or in anti-phase. The behavior persists 
even when the system is perturbed by noise (data not 
shown), evidencing the robustness of the phenomenon. 
In Fig. [U the inhibitor diffusion was varied in order to 
change the type of phase locking. As we show below, 
other parameters can be similarly tuned to control the 
system's behavior. 

Theoretical analysis in d = 1. The effect of cou- 
pling on the behavior of the excitable membrane can 
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Figure 1: Oscillations of the membrane inhibitor v in a 
one-dimensional bulk, bounded by two 0-dimensional "mem- 
branes". The left panel shows time series in each of the 
two boundaries (denoted respectively by a and b). Two 
regimes are shown: in-phase and anti-phase oscillations, la- 
beled "1" and "2", and are obtained by setting Dv = 0.7 
and Dv = 0.3, respectively. The two right panels show spa- 
tiotemporal representations of the inhibitor concentration in 
color code, with time running from top to bottom and space 
represented horizontally. The parameters used in these sim- 
ulations are L = 10, e = 0.015, k v = 5e, cry = Dv/W0 and 
era — Djj = k u = 0. As initial condition, the system is set at 
the rest state with additional (small) spatial fluctuations, so 
that there is a slight heterogeneity in the initial concentration 
field. 



be determined by using the solution and boundary 
conditions of the bulk equations (HJ-© m t ne mem- 
brane equations ©-(HI). To that end, we determine 
the stationary solutions of the diffusion equations ([T])- 
@ in a one-dimensional region of length L. In 
the case of the inhibitor, the resulting density pro- 



file is V(x) = V(0) cosh 



x) I cosh 



where V(0) = v. Its derivative at the boundary is then 



V'(x = 0) = -wj-^tanh 




(7) 



This expression, when inserted into the steady-state 
membrane equation 



= eg{u,v) + k v V'(x = 0), 

leads to a new effective inhibitor nulcline 

z u 

V = 



1 + b- tanh 

to 



L/2 



(8) 



(9) 



Here U = — and In = \/^ iL - Note that the initial in- 
hibitor nulcline v — zu is modified due to the coupling 
with the bulk, effectively decreasing its slope z as the cou- 
pling increases. The same analysis can be done for the 
activator, although in that case, the contribution from 
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Figure 2: Qualitative scheme of the effect of the bulk on the 
nulclines of the excitable membrane. A full circle denotes a 
stable fixed point and an empty circle an unstable one. The g 
nulcline labeled go corresponds to the absence of coupling with 
the bulk, gi and (72 define modified u-nulclines for increasing 
k v 7^ 0. The case of gi corresponds to a self-sustained oscilla- 
tory regime, while for larger coupling (52) the system is back 
at a steady state again. The / nulcline is practically unmod- 
ified by the coupling, and is left unchanged in this schematic 
plot. 



Figure 3: Phase diagram of a one-dimensional diffusive and 
degrading bulk delimited by two point excitable boundaries. 
Solid lines have been determined analytically. The oscilla- 
tory region lays within these lines. Symbols correspond to 
numerical simulations: full circles, stationary state; empty 
circles, in-phase oscillations; empty stars, anti-phase oscilla- 
tions; empty squares, independent oscillations. Parameters 
are those of Fig. [1] except for Dy = 0.5, and k v and L are 
varied. 



the coupling with the bulk barely modifies the shape of 
the membrane w-nulcline [f(u, v) = 0)], and therefore its 
effects will be ignored in what follows. 

With the above considerations, we can now under- 
stand, both qualitatively and quantitatively, the effect of 
membrane-bulk coupling. As shown in Fig. [21 the effect of 
the diffusive and degrading bulk can be mapped to an ef- 
fective variation of the slope of the g(u, v) nulcline, yield- 
ing a change from an excitable situation where only one 
stable fixed point exists, to an oscillatory regime where 
the fixed point is unstable and a limit cycle (not shown) 
develops. Finally, large enough coupling even leads the 
system back again to an excitable regime. The result ob- 
tained in Eq. ([9]) shows that three length scales control 
the system's behavior in the steady state: L, l and l\. 
The first one, L, is the natural length of the system. The 
second one, Zq, is a characteristic length determined by 
the ratio of the diffusion Dy and the degradation cry, 
and corresponds to an action length of the bulk. Finally, 
the effect of the coupling enters directly through the com- 
bination of the coupling strength k v and the time scale 
ratio e, giving rise to a characteristic length scale l±. 

Making use of Eq. © and considering the location and 
the stability of the fixed point (determined by the cross- 
ing of the nulclines), we can analytically predict the re- 
gion in parameter space where the system will behave as 
an oscillator. The corresponding phase diagram is shown 
in Fig. [H The theoretically predicted boundaries of the 
oscillatory region (solid lines in Fig. [3J agree very well 
with simulations. Therefore, the steady state analysis 
yields a clear understanding of why the excitable mem- 



brane becomes unstable and starts to oscillate. The sec- 
ond set of important properties of the system, i.e. syn- 
chronization and phase locking between the oscillating 
poles, have to be determined numerically. Such an inves- 
tigation leads to different types of oscillatory regimes, as 
shown in Fig. [31 in which anti-phase oscillations (stars) 
separate the region of in-phase oscillations (empty cir- 
cles) from a domain where the two poles of the system os- 
cillate independently (squares), namely, having the same 
frequency but an undefined phase relation. Note that for 
large L the pole oscillations can not interact, since the op- 
posing membranes are too far away. On the other hand, 
for large enough k v (i.e. li) oscillations (when they occur) 
can only be in phase. Finally, we note that increasing the 
inhibitor diffusion coefficient Dy makes Iq shorter, which 
changes the overall scale of the phase diagram. 

Extension to d = 2: numerical simulations. Once the 
one-dimensional case is understood and fully character- 
ized, we now consider a two-dimensional system with a 
rectangular geometry. The main difference is that now 
the membrane is also distributed in space (but with- 
out internal diffusion Q), while in the previous case it 
was zero-dimensional. However, by analogy to the one- 
dimensional case the relevance of the different parameters 
can be understood, and again the same phenomena are 
encountered. Every parameter has an important physi- 
cal meaning. Dy is responsible for the indirect coupling 
between membrane elements and is thus a necessary in- 
gredient for the synchronization. D\j is crucial to avoid 
several pathologies that occur in d > 1, such as accumu- 
lation of inhibitor along the membrane, which leads to 
local blocking of the oscillations, cry is key to the spon- 
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Figure 4: Temporal series of snapshots of V in 2- dimensions. 
A. Steady state solution. B and C. In-phase periodic oscilla- 
tions. D. Anti-phase periodic oscillations and traveling wave. 
The values of the parameters are L\ = 15.5, L1/L2 = 3.4, 
D v = D v = 0.1, a v = 0.00124, a v = 10""W, k u = k v = k 
and e = 0.015. The different regimes are obtained varying k: 
A, k = 0.03; B, k = 0.08; C, k = 0.095; and D, k = 0.11. 
Each snapshot is taken every At = 10. 



taneous onset of local oscillations at the membrane, as 
we have already seen in the one-dimensional case. k v de- 
termines the new effective nulcline due to the influence 
of the bulk, and controls the synchronization phases. 

Furthermore, there are now two length scales in the 
problem associated with the system size. L\ is the length 
of the rectangular system and its width is denoted by 
L2- Assuming L\ > L2, we are interested in observing 
oscillations along L\, mimicking the pole-to-pole oscilla- 
tions in elongated bacteria. In this case, the lateral walls 
can dominate the main activity of the system, preventing 
anti-phase oscillations for short lengths. For very long 
lengths the communication is fragile. For intermediate 
lengths Li, given that Li is smaller than L\, we expect 
the oscillations along the minor axis to be more or less 
synchronized in phase, and not to disturb excessively the 
pole-to-pole oscillations. 

We now present the phenomenology that can be ob- 
served numerically in two dimensions. In Fig. [51 tempo- 
ral sequences of two-dimensional snapshots of the density 
map of V are displayed for different kinds of behavior. 
The system is initiated again from the rest state, super- 
imposed with small random heterogeneous perturbations. 
When dynamical noise is added to the simulations, not 
only is synchronization maintained, but also the different 
regimes develop easier and faster. FigureUjA corresponds 
to small couplings, for which the system is unable to os- 
cillate and remains in a stable fixed point, corresponding 
to the quiescent state of the excitable membranes. For 
larger coupling (panel B), in-phase oscillations emerge 
from the center of the system. Panel C shows a differ- 
ent kind of in-phase oscillation, which is governed by the 



poles of the system. Finally, panel D shows a traveling 
wave that leads to anti-phase oscillations of the poles, 
similarly to what was found in the one-dimensional case. 
The traveling wave alternates periodically and slowly its 
direction of motion, from left to right and vice versa (with 
a period larger than the time span shown in Fig. [5J). In 
spite of these alternations, the poles oscillate periodically 
and in perfect anti-phase. The dimcnsionlcss parameters 
leading to this regime, given in the caption of Fig. [4j cor- 
respond to reasonable biological values when turned into 
dimensional units. For instance, just by assuming a bac- 
terium of length L ~ 4/im, it is found that Iq — 2.28fim 
and li — 1. 88/wn. 

Conclusions. We have reported a minimal mech- 
anism that generates pole-to-pole oscillations in non- 
active elongated media. The mechanism relies on the 
interaction of the system's bulk with an excitable (non- 
oscillatory) membrane, and can be understood analyti- 
cally for one-dimensional bulks making use of a phase- 
plane picture of the membrane's excitability. Both in- 
phase and anti-phase oscillations can be observed. In 
the case of a two-dimensional bulk, anti-phase dynamics 
is associated with a traveling wave that periodically re- 
verses its propagation direction. This could be a generic 
mechanism leading to spatiotemporal oscillations in sys- 
tems limited by active boundaries, such as cells. 
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